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The wave equation of low-frequency density waves in Bose-Einstein condensates at vanishing 
temperature in arbitrarily anisotropic harmonic traps is separable in elliptic coordinates, provided 
the condensate can be treated in the Thomas-Fermi approximation. We present a complete solution 
ON ■ of the mode functions, which are polynomials of finite order, and their eigenfrequencies which are 

characterized by three integer quantum numbers. 
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, 1 I. INTRODUCTION 

> ' 

Bose-condensates differ from normal gases or fluids by the existence of a macroscopic wave-function, their order 
parameter. The macroscopic wave-function deeply influences the spectrum of low-lying elementary excitations in 
Bosc-condensed systems, which become collisionless density waves with the velocity of sound. In the Bose-Einstein 
condensates of alkali-metal vapors in traps these collective sound waves have a discrete spectrum which is determined 
by the shapes of the trapping potential and of the condensate. Many experimental and theoretical JI-12 studies 
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^\ , have been devoted to their study. For Bose-Einstein condensates at zero temperature which are sufficiently large to 
validate the Thomas-Fermi approximation Stringari |j] found an analytical solution for the sound modes and their 
eigenfrequencies for the case of spherically symmetric parabolic traps. In the same work he even determined some 
of the eigenfrequencies for axially symmetric anisotropic traps. In subsequent works jl0|-|l^] the complete solution 
for the axially symmetric anisotropic case was given. In particular, it was demonstrated that the axially symmetric 
anisotropic problem forms a completely integrable system by exhibiting explicitly a third conserved operator B besides 
the wave-operator G and the axial angular momentum operator L z . The eikonal or 'classical limit' of the sound-waves, 
£j , determining their characteristic rays, was also studied in [ fl2| in the axially symmetric case. In this 'classical limit' 
I* ■ also the completely anisotropic case of a triaxial harmonic trap was investigated . It was shown that even in this 
. £h ! case the wave operator in eikonal approximation remains separable in elliptic coordinates. The complete integrability 
was demonstrated by exhibiting three phase-space functions G, B and A in involution. However, the solution for the 
mode- functions and eigenfrequencies was not yet given in |T^ ]. 

In the present paper we wish to return to this completely anisotropic case. ^From a practical point of view this has 
become of interest, because the first completely anisotropic trap has now appeared on the experimental scene p3[ . 
There the reported trap-frequency ratios are uj\ : uj\ : uj\ = 1 : 2 : 4. We shall return to this case when we give a 
numerical example at the end. 

^From a theoretical point of view the problem is also of considerable interest. The previous results on its classical 
limit suggest that also the full wave-operator remains separable in the general anisotropic case. This is indeed the 
case, as will be shown here. In fact we shall see that this problem can be related to a novel class of completely 
integrable elliptic billiards on an inhomogeneously curved space of arbitrary dimensionality. The sound modes in the 
Bose-Einstein condensate correspond to the 3-dimensional quantized version of such a billiard, their characteristics 
or rays are given by their classical limit. Furthermore it turns out that the classical limit of the billiards (in arbitrary 
dimension) can be connected mathematically to the equations of motion of an integrable system first studied by 
C. Neumann Q 140 years ago: a mass point on the sphere \x\ = 1 under the influence of an anisotropic harmonic 
force. 

The problem of collective modes in Bose-Einstein condensates can be viewed as the problem of small perturbations 
of the macroscopic wave function around its static equilibrium. At zero temperature the macroscopic wave function 
satisfies the Gross-Pitaevskii equation |T|] 
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where 



-^y 2 + U(x)- fi+^^l^x,^ 2 U(x,t) =ih<t>(x,t) ! I.I) 

U(x) = ^m(uj 2 1 x 2 1 +u J lxl+Ujlxl) (1.2) 

is the anisotropic harmonic potential of the trap, /i is the chemical potential, fixed by the requirement that N — 
J d 3 x\<fi(x, t)\ 2 , and ao is the s-wave scattering length, assumed to be positive throughout this paper. The macroscopic 
wave function (at zero temperature) is related to the number density n(x,t) — \<ft(x,t)\ 2 and the momentum-density 
g = ^r(0*V(^i — <j)V(f)*) — n(x)v s which satisfy the euqations of motion 

^ + V-nv s = (1.3) 
dv s r l 2 h 2 V 2 Vi 1 47^00 o 



We shall assume that the condition NaQ\/ muj/h 3> 1 is satisfied, where w = (wic^a^) 1 ' 3 , so that the Thomas- 
Fermi approximation can be applied to (1.3), where the term proportional to (V 2 i/n)/v / ^ is neglected. Then 



the equilibrium solution with dn/dt = = dv s /dt is given by v s = 0, n = no(x) = in $ a (/i — U(x)), /1 
(hu>/2) (lhNao/d) 2/5 , where J = 

tocj. Thus the condensate forms a triaxial ellipsoid. The collective excitations 



in the same approximation are now determined by eqs. fll.q ), linearized around the equilibrium solution 

n = uq + Sn, v s = 5v s 

V • n a (x)Sv s = (1.4) 



dt 

d5v s 4nh 2 ao 
dt m 2 



VSn = . 



Eliminating Sv s and with the ansatz Sn(t) — ip(x)e luJt we are left with the wave equation, first derived along the 
present lines by Stringari || 

^ = -V-(l-||-f-|w. (1.5) 
Here we introduce the characteristic lengths 



^, b =J^, c= K U», (1.6) 



which are the three semi-axes of the condensate ellipsoid. We note that 



2 1 1 1 

LUr, = — : — : — . (1.7) 

A a 2 b 2 c 2 K ' 



We also introduced the velocity of sound Co = wji/m in the center of the trap. In the following we assume a 2 > b 2 > c? 
without restriction of generality, i.e. u>\ is the smallest of the three trap frequencies and the largest. We sha ll 
sometimes use the notation a\ — a, a 2 = 6, a 3 = c. In this paper we shall be concerned with the solution of eq. fll.5|). 



II. EIKONAL APPROXIMATION 



To get the eikonal approximation to the wave equation (1.5) we first return to its explicitly time-dependent form 



replacing to 2 
Hamiltonian 



-d /dt 2 , and replace space and time derivatives via i 
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H 
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H, —ih"V —> p. This leaves us with the 



(2.1) 



whose trajectories describe the characteristics of the wave equation from which it was derived. Some features of the 
classical dynamics described by eq. (2.1) where studied in ref. p"2| ]. Here we wish to make a number of additional 
points. 
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a. Connection to Neumann's system 

The Hamiltonian equations of motion following from eq. (2.1) (with time now measured as a length by taking 
units with c = 1) 



1 - 



1 - 



■'"2 
b 2 



(2.2) 



v 2 



have the interesting property, remarked in ||12| , that the dynamics of the unit vector p = p/p can be decoupled 
from the dynamics of p. In fact, eliminating Xi from eqs. (2.2) we obtain 




and the equations of motion 



Pt = Pi 




(2.3) 



with the constraints p 2 — 1, p ■ p = 0. These are formally the equations of motion of a particle with unit 
mass with coordinates p on a sphere p 2 = 1 under the influence of the force F with components Fi = pi/a 2 
which were studied by Neumann p"4| . A discussion of this problem within the modern mathematical theory 
of integrable systems has been given by Moser jlT]]. He derives its conservation laws by constructing a matrix 
whose eigenvalues are preserved under the dynamics (|2.3|) . The conserved quantities (see also [^8|) are 



M k 



-Pk 



^ a 2 

i^k 1 



'—^ifiiPk ~ PkPif 



(2.4) 



They satisfy J2k^ k = — P 2 = We n °te that the Hamiltonian H (2T) is no longer amo ng these integrals 
because the abso lute value p of the momentum was eliminated in the derivation of (|2.3| ) from (2.2). However the 
conservation of (^J) can be used to recover the motion of p from the solution of ( p.3| ). There is a new obvious 
'energy'-integral of eq. (|2.3|) which is given by 



p 2 A l^M, 



It can be expressed in the original x, p variables as 

En = — 



1 \ - J P± ( , _ V^ffc^ _ A 2 , X iPi \ - X kPk \ 

2H 2 ^\a 2 \ l ^4) af + a 2 % a 2 ] 



(2.5) 



(2.6) 



i.e. it is now a quite complicated looking and far from obvious first integral of eqs. (2.2). With so me la bor it 
can be expressed in terms of the first integral A which will be introduced in sections 2c and 3 (see (3.16)) as 



E N = 



1 



.4 



3 1 

T 2 



(2.7) 



Another simple linear combination of the M% is 

B N = - ^ a l M k = XI a kPk + ^Y1 °* °* (PiPk ~ PkPi) ' 



It can be expressed in terms of the first integral B introduced in section 2c and 3 (see ( 3.1 6| ) ) as 



Bn — 



^{J2 a lpl + (vp)) 2 = 



H 2 



B 
H 2 
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b. Connection to a billiard on a curved space 

The Hamiltonian H describes the geodesic motion of a particle in a space with the metric 



ds 2 



1 - Ej -A 



This metric is inhomogeneous and conformal to the Euclidean metric. Its Riemann scalar curvature is given by 



^1 a?. 



and singular on the surface-ellipsoid. The distance to the surface-ellipsoid . x 2 /a 2 = 1 from any point inside 
the surface is finite. Furthermore, the metric velocity ds/dt is constant and simply given by ds/dt = 1 in 
our present units. The billiard particle always reaches the surface-ellipsoid perpendicularly and with diverg- 
ing orthogonal and finite tangential components of the momentum and is reflected with conserved tangential 
momentum component p"2| . 

c. Separation in ellipsoid coordinates 

Let us introduce elliptical coordinates E^] X,fi, v as the three roots of 



= 1 



a 2 + g b 2 + g c 2 + g 
We order these roots p = X,fJ,,v according to 

- a 2 < v < -b 2 < fi < -c? < X < . 



(2.9) 



(2.10) 



Geometrically surfaces A = const are ellipsoids, [i — const are one-sheeted hyperboloids, and v = const are 
two-sheeted hyperboloids, all confocal to the basic ellipsoid 

2 2 9 

— - A i- A i = 1 

a 2 b 2 c 2 

Explicitly, the Xi are given in terms of the new variables by 



Xl = ±4 



V + A)( a 2 +/i)(a 2 + ^ 



and cyclic. 



(a 2 - b 2 )(a 2 - c 2 ) 
The Hamiltonian-Jacobi equation 

uj = H{VS,x) 

is then separable, as shown in mM. Written in elliptical coordinates it takes the form 







(/i - v) 



(a 2 + A)(6 2 + A)(c 2 + A) 



dS_ 
dX 



a 2 b 2 c 2 u 2 



4A 



+ cyclic 



which is satisfied only if the angular bracket is equal to a linear function of A, i.e. 

2 -,2l2^,2, ,2' 



(a 2 + X)(b 2 + X)(c 2 + X) 



dS_ 

9A 



4A" 



1 



(—A — BX) and cyclic 



(2.11) 



(2.12) 



(2.13) 



(2.14) 



where A and B are separation constants. Eqs. ( [2.14 ) can be solved for A,B,u> 2 = H 2 as functions of the 
coordinates A, /z, v and the canonically conjugate momenta p\ — dS/dX, p M = dS/d^i, _£v — dS/dv. The 
results, expressed in terms of the Cartesian coordinates and momenta have been given in |12] and need not to 
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be given here. In any case they are easily recovered from the operators A, B derived in section III upon taking 
the classical limit — i?iV — > p, h — > 0. Here we wish to remark on a direct geometrical significance of the values 
of the three conserved quantities H 2 , A, and B. Let us introduce to this end the roots Ai > A 2 



of the quadratic equation 

BX 2 + A\ 2 + C = (2.16) 
where C = a 2 b 2 c 2 uj 2 > 0. The equations for the momenta p\, p^, p v now take the form 



/ B(X - Ai)(A - A 2 ) , , 

PX = ± \j- X(a 2 + X)(b 2 + X)(c 2 + X) andc y clic - ( 2 - 17 ) 

For a physical motion inside the ellipsoid each of the momenta p\, p^, p v needs to have an outer and an inner 



turning point. Thus there must be six turning points, which, according to eq. (2.17) can only be at values Ai, 
A 2 , 0, —a 2 , — b 2 , —c 2 . Of these the turning points at —a 2 , —b 2 , — c 2 correspond to coordinate singularities, 
the turning point at is the surface of the condensate and the turning points at Ai, A 2 correspond to caustic 
surfaces. Therefore Ai, A 2 must be real, and negative in order to be in the range of A, fi, i>, which imposes the 
conditions 

A 2 >ABC, A>0, B>0. (2.18) 

There are then four possible cases for the roots Ai, A 2 in which p\, p M , p v are real. These are: 

1) -a 2 < v < -b 2 < A 2 < fi < -c 2 < Ai < A < 

In this case p\ turns at p\ = on a surface of inner turning points forming the ellipsoid A = Ai and 
similarly p M turns at p^ — on a second surface of inner turning points forming a one-sheeted hyperboloid 
/i = A 2 . 

2) -a 2 < v < -b 2 < A 2 < fi < Ai < -c 2 < A < 

Here p^ turns at an outer surface fi = Ai and an inner surface fi = A 2 which are both one-sheeted 
hyperboloids. 

3) -a 2 < v < A 2 < -b 2 < fi< -c 2 < Ai < A < 

Here p\ turns on an inner ellipsoid A = Ai and p v turns on an outer two-sheeted hyperboloid v = A 2 . 

4) -a 2 < v < A 2 < -b 2 < fi < Ai < -c 2 < A < 

Here p M turns on an outer one-sheeted hyperboloid fi = Ai and p v turns on an outer two-sheeted hyper- 
boloid. 

The conservation of A and B for given to and the particular value taken by these quantities thus is reflected in 
the geometry of the two caustic surfaces occuring in each case. Similar results have been obtained for billiards 
in Euclidean space with ellipsoidal boundaries p(| . 

d. Semiclassical quantum numbers 

A 'semiclassical' approach to the soltion of the wave equation could be the application of the Bohr-Sommerfeld 
rule 

I\ = <L p\dX = n\ + 1/2 and cyclic. (2-19) 

27T ./ 



See ref. |21| for the case of isotropic traps. Here the quantum numbers n\, n^, n v can be interpreted to count the 
number of nodal surfaces with A = const, fi = const, v = const, respectively, which are ellipsoids, one-sheeted 
hyperboloids, and two-sheeted hyperboloids respectively. The original conserved quantities A, B and lo 2 can 
be expressed in terms of these quantum numbers and are thereby quantized in terms of the three independent 
integers n\, n lll n v . In section IV these quantum numbers will reappear in the exact solution of the wave 
equation in the slightly different notation 71,3, n 2 , n±, the correspondence being n\ — 71,3, ~ ?i 2 , n v = n\. 
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III. SEPARATION OF THE WAVE EQUATION 



Let us now return to the wave equation 

2 2 2 

cV = -v-(i-^-g-^w (3-1) 

a o c 

where we again adopt units with Co = 1. We are interested in the solutions ip of eq. ([O]) in the Hilbert space with 
the scalar product 



(^IVj) = / d 3 x^*(x)^(x) (3.2) 



'3 / * 

"' ' ' J 
TF 



where the integration is extended over the interior of the Thomas-Fermi ellipsoid. With this choice of the scalar 
product the operator G = — V • (1 — ^ x 2 /a 2 )V is self-adjoint if we pose as boundary condition 

on the boundary, i.e. the normal derivative d n ip should diverge there less than (1 — y]„. x 2 / a 2 ) -1 . This can be satisfied 
by choosing the Hilbert space of polynomials of finite order with the scalar product (f^J). These polynomials can be 
divided into eight different parity classes depending on whether they are even or odd under the inversion of any of 
the three coordinates x\, x%, X3. To be specific we put 



where a, (3, 7 = 0,1 determine the parity class and P m is a polynomial of order m in x\, x?,, x\. Eq. ( 3.1 
transformed to elliptical coordinates where it takes the form 

with 

F(Q) = (a 2 + g)(b 2 + q)(c 2 + e ) ■ (3.4) 



Eq. (3.3) can be solved by separation of variables via 

tp = <p\(\)<pn(fj,)(p v (v) (3.5) 

which leads to 

= 0*- u)g x (\) + (y- X)gM + (A - n)g v {v) (3.6) 

with 

, \ f r 7—] / \ d 2 <F(q) 1 . d a 2 b 2 c 2 uj 2 -. 



If eq. (ETa) is rewritten as 



(A)a j#i^ +A fe(!j (3 . 8) 

fj, — u n — v 



it becomes apparent that it can hold as a identity in A, /1, v only if 

Sa(A) = -4 - ?A 



4 4 

5^) = -j - (3-9) 



G 



where A and B are the same constants in all three of eqs. (3.9). ^,From eq. (3.£) and (3.7) the three constants us 2 , 
A, B can be expressed as eigenvalues of certain corresponding operators G, A, B. To do this explicitly we define the 
operator F e for arbitrary g as 



because then eqs. (3.9) with ( [3.7] ) can be rewritten simply as 

p g ip = {a 2 b 2 c 2 LU 2 + Ag + Bg 2 )iP 



(3.10) 



(3.11) 



where ip is the total wave function and g = A, /i, v. Solving for Lu 2 ip, Aip, Bip we find the simultaneous eigenvalue 
equations for u 2 , B, A, namely eq. Q which we abbreviate as lo 2 i)j = Gip and 



Bip = Bip. 



Aip = Atp 



with 



B = 
A = 



(A - - u)(v - A) 
4 



{(/* - v) [AF(A)^ + (F(A) + \\F\X)) ±] + cyclic^ 



(A - ,)(, - „)(„ - A) 1 ^ - ^ ^W 2 + + l XF '^ + C ^ C 



A lengthy but straight-forward calculation then yields the operators A and B in Cartesian coordinates 



A = | [(b 2 + c 2 )(x 2 - a 2 ) + a 2 (x 2 + x 2 )} ^ + 2a' 



if 







X2X3 



+ 3(6 2 + c 2 )xi-^— + cyclic 

0X2OX3 OX\ 



r) 2 B 2 
2 .2 ° 



B = (x.V) (a :.V + 3)-a^-^ . 



(3.12) 

(3.13) 
(3.14) 

(3.15) 
(3.16) 



By construction the eigenvalue equations for A, B and G can be satisfied simultaneously, i.e. these operators must 
commute, as one may also check by explicit calculation 

[G,A] = [B, A] = [G,B] = 0. 

In the axially symmetric case, e.g. b 2 — c 2 , the operator A may be expressed in terms of the angular momentum L z 
around the axis of symmetry, here chosen as the 1-axis, according to 



A = c l B + {a 2 - c 2 )L 2 z + a 2 c £ G 



,2„2/ 



IV. SOLUTION OF THE WAVE EQUATION 



After the separation of variables the equation to be solved follows from eq. (3/7) with (3.9) as 



[ - QF(e)j2 - (Hq) + 2 gF ' {g)) T g + i(« 26 w + a q + V)] Mq) = 



(4.1) 



for g = A, /1, v, i.e. precisely the same equation appears in all three elliptical coordinates. We now restrict the solutions 
of fl4.ip to the space of polynomials in Cartesian coordinates. This means that in elliptical coordinates they must be 
of the form 



Vq{0) = |a 2 + el f \b 2 + g\*\c 2 + g\iP m (p) 



(4.2) 



where P m ( p) is a polynomial of order m and the exponents a, 0, 7 can take o n th e values and 1. Inserting this ansatz 
in eq. (4.1) results in an expression containing the same prefactors as eq. (4.2) but multiplied with a polynomial of 
order m + 2 whose coefficients must all vanish, yielding m + 2 equations from which the three eigenvalues A, B, w 2 and 
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the m — 1 unknown coefficients of P m (one coefficient is fixed by normalization) must be determined. The vanishing 
of the coefficient of highest order fixes the value of B as 

B= (2m + a + /5 + 7 )(2m + a + /5 + 7 + 3). (4.3) 

The condition that the coefficient of the next to highest order term vanishes in principle fixes the value of A. However, 
it turns out that A, an d al so uj 2 , cannot be determined without at the same time calculating all the coefficients of the 
polynomial P m in eq. (4^2). Simple results are therefore not obtained in this way, except for a few of the lowest lying 
modes. 

Therefore a different procedure is used. It is a generalization of the analysis used in the solution of the Lame 
equation |]l9|. We shall assume now, and make plausible at the end of this section that the eigenvalues of the three 
separation constants uj 2 , A, B uniquely specify the corresponding eigenfunction, up to an arbitrary multiplicative 
factor. Then one can choose 

tp\(o) = <p»{q) = <Pv(o) ■ ( 4 - 4 ) 



In the polynomial ansatz (4.2) we write 

m 

Pm(Q)=H(Q-0i) (4.5) 



where Oi are the (possibly complex) roots of the polynomial P m . The ansatz for the total wave function ip 
ipx(^) i Pfj,(l^) i P^(i / ) then becomes 

ip = const\{a 2 + A) (a 2 + /i)(a 2 + is)\ a/2 \(b 2 + \){b 2 + y){b 2 + v)fl 2 



■2 



Using eqs. (2.11) and the identity 

„2 



A)(c 2 + M )(c 2 + v)\~<' 2 jj(A - e t ){n - 6i)[y - B t ) (4.6) 



Oi b 2 + 9i c 2 +9 t (a 2 + 9 t ){b 2 + 9A(c 2 + 0A 



the wave function (4.6) can be written rather simply in Cartesian coordinates as 

rn / 2 2 2 \ 



Oi ' b 2 
i=u - 

where for i = the factor under the product is defined as 1. In ( |4.8| ) and in the following we omitt a normalization 
constant and work with unnormalized wave functions. They are completely parametrized by the yet unknown pa- 
rameters Oi which dete rmine the nodal surfaces. As can be seen directly from (^?j) the remarkably simple form of 
the wave function ([D]) is a direct consequence of the separation ansatz. For all of the following considerations the 
form (4.8) of the wave function will be used. It follows from ( |4.8| ) that the nodal surfaces are quadrics confocal to 
the Thomas-Fermi ellipsoid. In order to relate the eigenvalues of G, A, B to the 6i we insert the ansatz ( |4.8| ) in the 
eigenvalue equations. The calculations become simpler with the use of the intermediate variables 

2 2 2 

= 4- 4- 3 _ i 

^ a 2 + 6i b 2 + 8, c 2 + 8 t 

2 2 9 

^o = ^ + | + ^f-l (4.9) 



a 



n=n^ 

i 

Then 



an^ =n d 2 n _ an 
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^From Gip = w 2- with i/j = x"x^ xJU we obtain 



*«f*i{g + * + J -4f;i]n 

k 4=1 



c»n , 4 4a + 2 4/3 + 2 4 7 + 2 w 



Thus the 9i have to satisfy 



4 t 4a + 2 ( 4/3 + 2 ( 4 7 + 2 

r2 I fl, + 2-^ 



i 2 + 6»i 6 2 + i 



= 



y 2 "J 



(4.11) 



for i = 1, 2, . . . , m. Here and in the following denotes the sum without the diagonal term. The eigenvalues of G 
become 



u? 2a 2/3 2 7 „ ^ 1 
— = 1 — — H — 4 > — 



(4.12) 



(where we momentarily restored c\ for later convenience). For m = the sum X)I=i nas ^° ^ e interpreted as 0. 
Similarly we obtain for B 



an 



= x^x^xjl (2m + a + /3 + 7 )(2m + a + /3 + 7 + 3)n + ^ — 6^(0) 
I i=i ^ 



(4.13) 



which leads again to the condition Gi(8) = and otherwise gives back the eigenvalue ( |4.3| ) for B, which, in particular 
turns out to be independent of the 0j. Finally applying A to ip we obtain 



Aip = xfx^xll [2/3 7 a 2 + 2 7 a& 2 + 2a/3c 2 



+ (4m + 3)(a(6 2 + c 2 ) + /3(c 2 + a 2 ) + 7 (a 2 + 6 2 )) 
+4(a 2 + b 2 + c 2 )m(m + 1) 

m 

+ (4a + 4/3 + 4 7 + 8m + 2) ^ 

i=l 

+E ^ t(« 2 + fe2 + c2 - »? - - 2 - + 02 ] Gi(9)\ . 



Again Gi{9) = must be satisfied, and the eigenvalues of A then can be read off eq. ( 4.14 ) in the form 
A = 4 (a + (3 + 7 + 2m + -) ^ (9; 

2 

+ [2a/3c 2 + (4m + 3)a(6 2 + c 2 ) + 4m(m + l)a 2 + cyclic}. 



(4.14) 



(4.15) 



It remains to determine the 9i from the m equations (4.11) Gi(0) = 0. First we show that all 9i are real. To prove 
this let us suppose they are complex, in which case they also satisfy Gi(9*) = 0. Hence 



o = '£(e i -e*)(G i (9)-G i (e*)) 



(4.16) 



which, after some algebra, and writing 0i = \9i\e *, leads to 



i 



4a + 2 4/3 + 2 4 7 + 2 2 . 2 

+ TTo ! 7TT7 + in.* ,o )|fli| sin 1?j 



|0 4 | 2 |a 2 + ^| 2 |6 2 + 4 | 2 |c 2 + ^| 2 ' 
(|^| sin^- 1^1 sin i?,) 2 



(4.17) 
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This condition is only satisfied if all phases #j equal either or ir, i.e. if the 9i are real. A permutation of two 9i does 
not lead to a new eigenfunction. Hence the 9i can be assumed ordered according to 



m < 9 m -i < ■ ■ ■ < 0i ■ (4.18) 
Let us consider the cases m = 0, m = 1. For to = we obtain from (4.12) eight eigenfrequencies 

for a,P,j= 0, 1. These are the frequencies lo 2 = 0, u)\,u)\, w|, o;f +a;|, o;| +W3, For m = 1 one 
has to solve just a single cubic equation 

4 4a + 2 4/3 + 2 4 7 + 2 s 

01 a 2 + 61 b 2 + 6 X c 2 + 0i v ' 

to find three different values for X , which we order according to 

-a 2 < 0i,3 < -6 2 < 0i,2 < -c 2 < 0i,i < (4.21) 

which correspond, for each of the eight choices of the tripel (a, j3, 7), to three independent solutions. Taken together 
the case m = 1 therefore gives 24 different frequencies 

4c 2 

lu 2 = awl + + juj 2 - — ^ jfe = l,2,3. (4.22) 

01, k 



For the special case a,/3, 7 all vanishing the result ( 4.22 ) with the cubic equation ( 4.20 ) was already obtained by 
Stringari p2^ |. The three roots 01^ correspond to wave functions with a single nodal surface, which for k — 3 is a 
two-sheeted hyperboloid, for k — 2 a one-sheeted hyperboloid and for k = 1 is an ellipsoid. 

Let us now turn to the case of general order m of the polynomial mode function. It is quite remarkable that in this 
case the m equations Gj(0) = can be derived as the extrema of a single potential V(9) 

G i (0) = -^M = O (4.23) 

where 

1 rn 1 m /Q 1 m 

8V(0) = -2$>N-(f + 3 )5>|a 2 + i |-(f + i )5>|6 2 + i | 

2= 1 z i 

1 rn mm 

i i—1 j=i-\-l 

Thus the problem has now become the following exercise in statics: In a one-dimensional space one has four fixed 
positive fictitious point-charges aligned along the negative 0-axis, namely 

a charge 1/2 at = 0, 

a charge 2- + 1/4 at 9 = — a 2 , 

a charge | + 1/4 at 9 = -b 2 , 

a charge % + l/4 at 9 = — c 2 , 

between these fixed charges to movable positive point-charges of unit strength, and interacting among themselves 
and with the fixe d cha rges with 1-dimensional inverse-distance forces, are t o be distributed in such a way that a 
force-equilibrium ( 4.23 ) may result. It is clear from the form of the potential ( 4.24 ) that the movable charges can be 
distributed arbitrarily over the three intervals, namely 

77,3 charges in — c 2 < 9 < 0, 

n 2 charges in — b 2 < 9 < — c 2 , 

ri\ charges in —a 2 < 9 < —b 2 

with rix + ri2 + «3 = to. There are ( m ^ 2 ) = "m\2\ wavs to make this distribution, each leading to a different unique 
equilibrium configuration for the movable charges at positions 9i . It is immediately clear from the mechanical analogy 
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that all 9i must be different. The three integer numbers m, n 2 , ^3 are therefore the natural quantum numbers of the 
problem. As the eigenvalue B is independent of the positions of the movable charges, it must be ( m ^~ 2 )-fold degenerate 
for each of the eight choices of the triple (a, (3, 7). The uniqueness of the equilibrium distribution of the 0, for given 
711, n 2 , ri3, which determine uniquely the eigenvalues lu 2 , A, B, justifies a posteriori the uniqueness assumption made 



after eq. (4.3). It is very easy to find the minima Oi of the p oten tial (4.24) numerically for any desired tri ple o f integers 



n\, n 2 , ri3 and to determine thereby the mode frequencies ( 4.12 ) and the corresponding mode functions (4.8). In table 
1 we give the 20 lowest lying mode frequencies for the experimentally realized case lo\ : uj 2 : uj\ — 1 : 2 : 4 in units 
of the smallest trap frequency lo\. Also given there are the parities a, /3, 7 (e.g. there is even parity in x\ if a = 
and odd if a — 1), and the quantum numbers 123, 77,2, rii- The latter give, respectively, the numbers of nodal surfaces 
of ellipsoidal, one-sheeted hyperbolic and two-sheeted hyperbolic form, all confocal to the Thomas-Fermi ellipsoid. 
The ellipsoidal nodal surfaces counted by are more elongated ellipsoids inside the Thomas-Fermi ellipsoid and 
would correspond to radial waves in the spherically symmetric case. As the 9- values counted by 713 are in the group 
-c 2 < 9 < with the smallest absolute values, it is clear from ( 4.1 2| ) that the ellipsoidal nodal surfaces lead to the 



highest frequencies. The one-sheeted hyperbolic nodal surfaces counted by n 2 are ellipsoidal hyperboloids around the 
X3-axis which intersect planes X3 = const in ellipses and planes x 2 = const • xi in the two branches of hyperbolas, the 
hyperbolas opening up in all directions orthogonal to the a^-axis. Finally, the two-sheeted hyperbolic nodal surfaces 
counted by ni are formed by the two branches of ellipsoidal hyperboloids around the 1-axis opening up in the positive 
and negative a:i-direction. They are cut in ellipses by planes xi — const and in hyperbolas by planes X2 = const and 
X3 = const. The quantum numbers n 2 and n\ are clearly the analogues of angular momentum quantum numbers (i.e. 
the quantum numbers of spherical harmonics) for the elliptic geometry. 

V. SYMMETRIC TRAPS AS LIMITING CASES 

The cases of axially symmetric and isotropic traps must of course be contained in our results as limiting cases. Let 
us see how. 

a. Axially-symmetric trap 

Let us suppose we have axial symmetry of the trap around the xi-axis. In this case we should study the limit 
b +c from above. This limit has to be taken in the expr ession for the wave function (|4.8|), in the force equation 



(4.11) and in the result for the mode frequencies ( 4.12| ). The positions Oi of the n 2 charges which have been 



distributed in the interval — b 2 < 9 < —c 2 all approach — c 2 in the limit. Therefore the mode frequencies in the 
limit become 



2a 2(3 + 2 7 + 4n 2 ^ 



i— ni+n2 + ly 



The force-equilibrium for the ri\ + n 2 charges outside the interval [— 6 2 ,— c 2 ], which alone enter the sums in 
eq. (EO), becomes 



4 4a + 2 4/3 + 47 + 8n 2 + 4 
fi + a 2 + 9, + c 2 + 9, 




(5.2) 



for i = 1, 2, . . . ,nx] n\ + n% + 1, . . . , m. The solution of eq. ( |5.2| ) is sufficient to determine the mode frequencies. 
However, the wave functions still depend on the asymptotics of the distribution of the n 2 charges in — b 2 < 9 < 
—c 2 . The force equations for i = m + 1, n% + 2, . . . , n\ + n 2 have of course to be handled with some care, as the 
interaction terms between these charges, which approach each other arbitrarily closely, diverge. To isolate and 
divide out the diverging prefactor, which would be automatically cancelled if we worked with normalized wave 
functions, we define parameters ti with — 1 < i; < by ni +i = —(? + ti(b 2 — c 2 ). Then the limiting term of the 
force equations in question for b — > c reads 



47 + 2 4(3 + 2 



n-2 
X — ^ 



0= — + — + , = 1,2,. ..,, 2 . (5.3) 

These e qua tions can now be solved by similar techniques as used before. Fortunately, however, it will be sufficient 
to use ( |5.3| ) without explicit knowledge of its solutions. Turning now to the wave functions we define cylinder 
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coordinates around the rri-axis by x\ = z, x 2 = gsunp, X3 = g cos (p. Using the variables ti and multiplying by 
a factor (b 2 — c 2 )™ 2 to eliminate the divergence the dominant term of the wave function for b — ► c becomes 



a? + 6i b 2 



(5.4) 



Hi m 



z a g 2 " 2 +/3+7 sin* 3 y> cos 7 <p Y[ (cos 2 V + ) II II 



z=l i=m— n3 + l 



c 2 + 0j a 2 + 1 



Using the force equilibrium (5.3) for the charges pinned between —b and — c it can be shown that 



sin* 3 9? cos 7 y> (cos 2 ip + U) ~ cos 



i=l 



(2n 2 +/? + 7 )<p-^ 



(5.5) 



As a result the wave function (5.4) for b — > c goes over to 



r/> - (? : 



(2n 2 + /? + 7 V- 



7T/3 



(5.6) 



ni m 

n n 



C 2 + I 



1 



2 — 1 2— m — 71-3+1 

The quantum number £ 2 for the conserved angular momentum around the z-axis can be read off from eq. fl5.6| ) 

|4| = 2n 2 +(3 + 1 . 



Clearly one may take linear combinations of the wave functions ( |5.6| ) for \£ z \ fixed but suitably changing /?, 7 
or n,2 to form eigenstates ~ e ±l '^' v of i 2 = —id/dip. 

b. Isotropic traps: 

We can finally take the further limit c — > a in the results of the preceding section. By similar considerations as 
described there we obtain for the mode frequencies 



J 2 2a + 2/3 + 27 + Ani + An 2 



m 1 

-4 £ * 



(5.7) 



u i=m— 713+1 

The positions of the ri3 free charges satisfy the force equilibrium 

4 4a + 4/3 + 4 7 + 8n a + 8n 2 + 6 , 







(5.8) 



E' 8 
6»i - fa 

3— m— TI3 + I J 



for i = m — n3 + l,...,m. 



Introducing g = rsini?, 2 = rcosi? the wave functions take the form 



V>~^p/ a| (costf)e^ J] 



- 1 



where 

£ = a + P + ~/ + 2n x + 2n 2 = a + /3 + 7 + 2(m - n 3 ) 
is the total angular momentum quantum number. Clearly 7^3 now is the radial quantum number. 



(5.9) 



(5.10) 
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Using the force equilibrium (5.8) it can be shown that 



n 



i— m— n3+l 

and, furthermore, that 



m 1 



i—m—n^-\-l 



m 1 



= n 3 n 3 



-m— 713 + I 



This simplifies the result (5.7) for the mode frequencies which now becomes 



uj 2 = lu 2 [27J3 + 2n 3 £ + 3n 3 + l] 
and is indeed the result originally derived by Stringari H , of course in a much more direct way. 



(5-11) 



(5-12) 



(5.13) 



VI. CONCLUSION 



We have provided a complete solution for the eigenmodes and eigenfrequencies of the anisotropic wave equation 
governing the low-frequency collisionless density waves in Bose-Einstein condensates in harm onic oscillator traps with 
arbitrary ani sotr opy in the Thomas- Fermi limit. The eigenfrequenc ies ar e given by eq.( 4.12| ), the mode functions are 
given by eq.( |4.8| ), where the parameters 9i are the solutions of eqs.(4.23), (4.24). The solution was possible, because 
the system was found to be completely integrable, with three mutually commuting operators G,A,B. However, 
unlike in many more familiar examples in quantum mechanics, the solution was not constructed by directly solving 
the simultaneous eigenvalue equations for the three commuting operators, because their eigenvalues, apart from that 
of B, turned out to be rather complicated, not providing us with simple quantum numbers. Rather our solution 
proceeded by first constructing a simple form of the total wave function which followed from the separation ansatz. 
The solution to the equations fixing the free parameters 9i in the wave function then provided us with the natural 
simple integer quantum numbers ni,?i2,n 3 of the problem, on which the eigenvalues of A and the mode frequencies 
depen d in directly and in a complicated way. We never even had to determine the eigenvalues of A explicitely. The 
form (4.8) of the wave function is remarkably reminiscent of the Bethe-ansatz. It might be interesting to follow this 
connection further as it might shed more light on the mathematical structure of the problem we have solved in this 
paper, and might e.g. allow us to gain a deeper understanding of the degeneracy of the operator B. The method 
of solution, and in fact even the detailed structure of the solution, generalizes directly to the analogous problem in 
an arbitrary number of dimensions.. Thus the physical problem we have considered is found to be a member of a 
whole family of integrable systems with connections, as we have discussed, to billards on a curved space conformal to 
Euclidean space, and the class of integrable systems discussed in the memoir of Moser p7[. 
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TABLE I. Mode frequencies for the trap J{ 


2 2 -i 

lj 2 : w 3 = 1 


2 : 4 in units of ui\. 


a, [3, 7 are the parity quant 


um numbers, 



ni,n2,ri3 are the three positive integer quantum numbers which label uniquely each mode function (see text). An accidental 
degeneracy occurs for the states (0, 0, 0, 0, 1, 0) and (1, 0, 1, 0, 0, 0) where uj/u>i — \/E. 
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